In [1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
plt.rcParams['figure.dpi'] = 200
plt.rcParams['figure.figsize'] = (4,3)

You might have heard a lot of statistics jargon surrounding PCA, such as "dimensionality reduction" and "maximizing variance". To build up to the meat of PCA, I want you to get used to thinking in terms of variance, vectors, and variance along those vectors. We'll start simple and build up from there.

First, consider the simplest case, of a one-dimensional vector $\vec{x}$, whose elements are just random numbers about some some mean value:

$$ \vec{x} = \mathcal{N}(\mu, \sigma) $$
In [2]:
x_dat = np.random.normal(5, 1, 10)

fig, ax = plt.subplots()
ax.plot(x_dat,'ko')
ax.set_xlabel('$i$')
ax.set_ylabel('x$_i$')
ax.set_ylim(0, 10)
plt.show()

There is only one dimension, $x$, along which the values $x_i$ can vary, so it's pretty easy to describe the variance of our sample here with only one number along the $x$ dimension, the mean $\mu$:

$$ Var(x) = \sigma^2_x = \frac{1}{N} \sum^N_{i=1} (x_i -\mu)^2 $$

Think about the variance as a function of $\mu$, and the sample mean $\bar{x}$ is the value of $\mu$ that minimizes the variance i.e.,

$$ \frac{d\sigma^2}{d\mu}(\mu = \bar{x}) = 0 $$
In [3]:
test_mus = np.linspace(0, 10)
test_vars = np.array([np.sum(  (x_dat - mu)**2) / 10  for mu  in test_mus])

plt.plot(test_mus, test_vars)
plt.xlabel('$\mu$')
plt.ylabel('Variance')
plt.show()

We can solve for this analytically by differentiating the expression for variance:

$$ \frac{d\sigma^2}{d\mu} = \frac{1}{N} \sum^N_{i=1} -2(x_i - \mu) = \frac{-2}{N} \Big[ \sum^N_{i=1} x_i - \sum^N_{i=1} \mu \Big] $$

Set that to 0 and solve for $\mu$:

$$ 0 = \frac{-2}{N} \Big[ \sum^N_{i=1} x_i - \sum^N_{i=1} \mu \Big] \rightarrow 0 = \sum^N_{i=1} x_i - \sum^N_{i=1} \mu \rightarrow \sum^N_{i=1} x_i = \sum^N_{i=1} \mu = N \mu \rightarrow \boxed{\bar{x} = \frac{1}{N} \sum^N_{i=1} x_i} $$

Hopefully this formula is familiar as the definition of the sample mean!

In [4]:
sample_mean = x_dat.sum() / 10

fig, ax = plt.subplots()
ax.plot(x_dat,'ko')
ax.axhline(sample_mean, color='r')
ax.set_xlabel('$i$')
ax.set_ylabel('x$_i$')
ax.set_ylim(0, 10)
plt.show()

Now let's consider some vector $\vec{y}$ of observed data as a function of $x$:

In [5]:
x = np.arange(-10, 11, 2)
m = 0.3
b = 3
y_dat = m * x + b + np.random.normal(0, 1, len(x))

plt.plot(x, y_dat, 'ko')
plt.xlabel('x')
plt.ylabel('y')
Out[5]:
Text(0, 0.5, 'y')

Similar to before, the data only varies along one axis - this time $y$ - but we can do better than just the mean to describe the data. The data here clearly looks like a typical line of the form $y = m x + b$. $x$ is a dependent variable here and is known. How do we solve for the values of $m$ and $b$? Just like the mean, we can solve for the values that minimize the variance:

$$ \sigma_y^2 = \frac{1}{N} \sum^N_{i=1} (y_i - (m x_i + b))^2 $$
In [9]:
test_ms = np.linspace(-1,1)
test_vars = np.array( [np.var( y_dat - (m * x + b_opt)) for m in test_ms])

plt.plot(test_ms, test_vars)
plt.xlabel('m')
plt.ylabel('Variance')
plt.show()

Let's solve for the optimal values of $b$ and $m$ analytically. First, let's differentiate with respect to $b$:

$$ \frac{d \sigma_y^2}{db} = \frac{-2}{N} \sum^N_{i=1} y_i - mx - b $$$$ 0 = \sum^N_{i=1} y_i - m x_i - b \rightarrow N b = \sum^N_{i=1} y_i - \sum^N_{i=1} m x_i \rightarrow b = \bar{y} - m \bar{x} $$

Now with respect to $m$:

$$ \frac{d \sigma_y^2}{dm} = \frac{-2}{N} \sum^N_{i=1} x_i (y_i - mx - b) $$$$ 0 = \sum^N_{i=1} x_i (y_i - mx - b) \rightarrow m \sum^N_{i=1} x_i^2 + b \sum^N_{i=1} x_i = \sum^N_{i=1} x_i y_i $$

We can plug in our previous expression for $b$:

$$ b = \frac{1}{N} \sum^N_{i=1} y_i - \frac{1}{N} \sum^N_{i=1} m x_i $$

so

$$ m \sum^N_{i=1} x_i^2 + \frac{1}{N} \sum^N_{i=1} x_i y_i - \frac{m}{N} \sum^N_{i=1} x_i^2 = \sum^N_{i=1} x_i y_i $$$$ N m \sum^N_{i=1} x_i^2 + \sum^N_{i=1} x_i y_i - m \sum^N_{i=1} x_i^2 = N \sum^N_{i=1} x_i y_i $$$$ N m \sum^N_{i=1} x_i^2 - m \sum^N_{i=1} x_i^2 = (N- 1) \sum^N_{i=1} x_i y_i $$$$ (N-1) m \sum^N_{i=1} x_i^2 = (N- 1) \sum^N_{i=1} x_i y_i $$$$ m = \frac{\sum^N_{i=1} x_i y_i}{\sum^N_{i=1} x_i^2} $$

Let's see how well that works:

In [7]:
m_opt = np.sum(x * y_dat) / np.sum(x**2)
b_opt = np.mean(y_dat) - m * np.mean(x)

plt.plot(x, y_dat, 'ko', label='data')
plt.plot(x, m_opt *x + b_opt,label='fit', color='r')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

Not bad!

Suppose we wanted to avoid all of this arithmetic, or the data was a bit more complicated and we wanted more coefficients to fit for, we can set this up as a system of linear equations:

$$ m \sum^N_{i=1} x_i + N b = \sum^N_{i=1} y_i $$$$ m \sum^N_{i=1} x_i^2 + b \sum^N_{i=1} x_i = \sum^N_{i=1} x_i y_i $$

Which can be expressed as:

$$ \begin{bmatrix} \sum^N_{i=1} x_i \ \ N \\\ \sum^N_{i=1} x_i^2 \ \ \sum^N_{i=1} x_i \end{bmatrix} \begin{bmatrix} m \\\ b \end{bmatrix} = \begin{bmatrix} \sum^N_{i=1} y_i \\\ \sum^N_{i=1} x_i y_i \end{bmatrix} $$

For reasons that will make sense in a second, let's call the first matrix $\mathbf{X}$, the second $\mathbf{B}$, and the third $\mathbf{A}$.

$$ \mathbf{X} \mathbf{B} = \mathbf{A} $$

We can solve for $\mathbf{B}$ by multipliying both sides by $\mathbf{X}^{-1}$:

$$ \mathbf{B} = \mathbf{X}^{-1} \mathbf{A} $$
In [10]:
X = np.array([[x.sum(),         len(x)    ],
              [ (x**2).sum(), x.sum()     ]])

A = np.array([ y_dat.sum(), (x*y_dat).sum()]).T

m_opt, b_opt = np.linalg.solve(X, A)

plt.plot(x, y_dat, 'ko', label='data')
plt.plot(x, m_opt *x + b_opt,label='fit', color='r')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

Even though just going through the arithmetic worked out fine, I wanted to introduce the idea that we can use various matrix operations to solve systems of equations and fit the data.

In general, we can assume that we can approximate most datasets $\mathbf{A}$ as a linear combination of n vectors $\vec{x}_i$:

$$ \mathbf{A} = \beta_0 + \beta_1 \vec{x}_1 + \beta_2 \vec{x}_2 + ... + \beta_n \vec{x}_n $$

where $\beta_i$ is the coefficient. A fairly common implementation is a polynomial fit, where each vector $\vec{x}_i$ is the dependent variable to the $i$th power:

$$ \mathbf{A} = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 + ... $$

But the vectors $\vec{x}_i$ can really be anything you want, like $\cos x$, a Gaussian, or really any kind of line that is appropriate for your data.

More succintly in matrix notation, we can put all the vectors we care about in a design matrix $\mathbf{X}$ where the columns are the vectors $\vec{x}_i$, and the coefficients in their own vector $\vec{\beta}$:

$$ \mathbf{A} = \vec{\beta} \mathbf{X} $$

In the case above with a simple line, the first vector will just be 1's and correspond to the intercept $b$, and the second vector is just the values of our dependent variable $\vec{x} = \begin{bmatrix} x_1, x_2, ..., x_N \end{bmatrix}$:

$$ \mathbf{X} = \begin{bmatrix} 1 \ \ x_1 \\ 1 \ \ x_2 \\ ... \\ 1 \ \ x_N \end{bmatrix} $$

The vector $\vec{\beta}$ just has two elements, the associated coefficients analogous to $b$ and $m$:

$$ \vec{\beta} = \begin{bmatrix} \beta_0 \ \ \beta_1 \end{bmatrix} $$

And then $\mathbf{A}$ in this case is just a single column vector of the observed $y$ values:

$$ \mathbf{A} = \begin{bmatrix} y_1 \\ y_2 \\ ... \\ y_N \end{bmatrix} $$

We then solve this linear system by taking the pseudo-inverse of $\mathbf{X}$ (we can't take the true inverse because it's no longer a square matrix):

$$ \boxed{ \hat{\beta} = (\mathbf{X}^\intercal \mathbf{X})^{-1} \mathbf{X}^\intercal \mathbf{A} } $$

This is called the Normal Equation and allows us to perform linear regressions. With our chosen design matrix and optimal value of $\beta$, we can fit the data pretty well within some random noise/error $\epsilon$:

$$ \mathbf{A} = \hat{\beta} \mathbf{X} + \epsilon $$
In [11]:
N = len(y_dat)
X = np.column_stack([np.ones(N), x]) #our design matrix

A = y_dat.reshape((N, 1)) #observed data, reshape into a column vector

B_opt = np.linalg.inv( np.dot(X.T, X)  ).dot( np.dot(X.T, A)) #normal equation

A_opt = X @ B_opt

plt.plot(x, y_dat, 'ko', label='data')
plt.plot(x, A_opt,label='fit', color='r')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

Recall that for matrix multiplication $\mathbf{A} \mathbf{B} = \mathbf{C}$, element $c_{ij}$ is the dot product (sum of the element-wise products) of the $i$th row of $\mathbf{A}$ and $j$th column of $\mathbf{B}$:

$$ c_{ij} = \sum_{k=1} a_{ik} b_{kj} $$

You'll find then with the Normal Equation and 1D data that the coefficient for any vector is the same for the expression for $m$ in the simple linear fit above:

$$ \beta_i = \frac{\sum_{j=1}^N x_{ij} a_{j}}{\sum_{j=1}^N x_{ij}^2} $$

Even for the intercept vector, this works out to the mean:

$$ \beta_0 = \frac{\sum^N_{j=1} a_j}{\sum^N 1} = \frac{1}{N} \sum^N_{j=1} a_j $$

This is the main idea behind PCA - any matrix $\mathbf{A}$ can be approximated as a linear combination of vectors $\mathbf{X}$. The question is then - how do we acquire the right vectors?

In the example above, $x$ was an exactly known dependent variable, but what if both $\vec{x}$ and $\vec{y}$ are randomly sampled variables?

In [14]:
mean = np.array([ 0,0 ])

cov=  np.array([ [1,0.55],
                 [0.55,0.5 ]   ])

x,y = np.random.multivariate_normal(mean, cov, 501).T

plt.plot(x, y, 'ko', alpha=0.8, markersize=2)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

There is a clear colinearity between $\vec{x}$ and $\vec{y}$ here, so we might be tempted to do a simple linear fit like above:

In [15]:
N = len(x)
X = np.column_stack([np.ones(N), x])

A = y.reshape((N, 1))

B_opt = np.linalg.inv( np.dot(X.T, X)  ).dot( np.dot(X.T, A))

A_opt = X @ B_opt

plt.plot(x, y, 'ko', label='data', markersize=2)
plt.plot(x, A_opt,label='fit', color='r')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()


plt.plot(x, y-A_opt.reshape(N), 'ko', label='data',markersize=2)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Residuals')
plt.legend()
plt.grid()
plt.show()

Well that's not very satisfying - there's still tons of variability in the data. However, would you believe me if I said this data set could be described by a linear combination of only two vectors, just like the simple line fit before? Enter Principal Component Analysis.

Part of our problem is that we are still thinking in terms of the data only varying along $\vec{x}$ or $\vec{y}$. But there are some latent/hidden vectors along which the data vary the most. These are the principal components. So rather than $\vec{x}$ or $\vec{y}$, we need to find vectors that $\textit{maximize}$ the variance, i.e., contain the most information.

Wait, we've been $\textit{minimizing}$ the variance up to this point, why do we want to maximize it now? Think of the above data set as $x-$ and $y-$ positions of people standing in a room, and you want to take a group picture. You want to maximize the amount of information about each person (their faces) so you want to come up with the vector of maximum variance and set up your camera perpendicular to that.

Let's test along a bunch of different angles for our new vector:

In [17]:
theta = np.linspace(0,  np.pi, 101) #angle from the x axis unit vector to rotate from
test_vectors = np.column_stack([ np.cos(theta), np.sin(theta)])
a = np.column_stack([ x, y]) #our data matrix, each column is a variable

rejections = np.array([  np.var(a.dot(b.T)/np.linalg.norm(b)) for b in test_vectors]) #calculate the variance based on vector rejections

plt.plot(theta, rejections)
plt.xlabel(r'$\theta$ [rad]')
plt.ylabel('Variance')
plt.show()

wheremax = (rejections == rejections.max()) #which angle rotation led to the highest variance?
#lets plot the new axis
b_max = test_vectors[wheremax][0]

plt.plot(x, y, 'ko', alpha=0.8, markersize=2)
plt.grid()
plt.plot([-b_max[0]*3, b_max[0]*3], [-b_max[1]*3, b_max[1]*3],color='b')
plt.xlabel('x')
plt.ylabel('y')
Out[17]:
Text(0, 0.5, 'y')

The way to do this exactly is to get the eigenvectors of the covariance matrix:

$$ Cov(x, y) = \sigma^2_{xy} = \frac{1}{N} \sum_{i=1}^N (x_i - \bar{x})(y_i - \bar{y}) $$$$ \mathbf{\Sigma}_{xy} = \begin{bmatrix} \sigma^2_x \ \ \sigma_{xy}^2 \\ \sigma_{xy}^2 \ \ \sigma_y^2 \end{bmatrix} $$
In [18]:
meanx = x.sum() / N #mean of x
varx = np.sum(( x - meanx)**2) / N #variance of x

meany = y.sum() / N # mean of y
vary = np.sum( (y - meany)**2) / N #variance of y

cov = np.sum( (x - meanx) * (y - meany) ) / N #covariance of x and y

#covariance matrix
K = np.array([[varx, cov],
             [cov, vary]])

val, vec = np.linalg.eig(K) #eigendecomposition, val are eigenvalues, vec is eigenvectors

origin = np.array([[meanx,meanx],
                   [meany,meany]])

plt.plot(x, y, 'ko', alpha=0.8, markersize=2)
plt.grid()
plt.quiver(*origin, vec[0], vec[1], scale=4, color=['b','r'], zorder=3)
plt.xlabel('x')
plt.ylabel('y')
Out[18]:
Text(0, 0.5, 'y')

As a sanity check, consider the limiting case of the very first example where we just had data varying along $\vec{x}$. Let's pretend it was actually a 2D data set, but the $y$ value is the same (0):

In [19]:
y0 = np.zeros_like(x_dat)

plt.plot(x_dat, y0,'ko')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

meanx = x_dat.sum() / 10
varx = np.sum(( x_dat - meanx)**2) / 10

meany = y0.sum() / 10
vary = np.sum( (y0 - meany)**2) / 10

cov = np.sum( (x_dat - meanx) * (y0 - meany) ) / 10

K = np.array([[varx, cov],
             [cov, vary]])

val, vec_simple = np.linalg.eig(K)

print('eigenvalues:',val)

origin = np.array([[meanx,meanx],
                   [meany,meany]])

plt.plot(x_dat, y0, 'ko', alpha=0.8)
plt.grid()
plt.quiver(*origin, vec_simple[0]*val[0], vec_simple[1]*val[1], scale=4, color=['b','r'], zorder=3)
plt.xlabel('x')
plt.ylabel('y')

eigenvalues: [1.10943516 0.        ]
Out[19]:
Text(0, 0.5, 'y')

The first principal component is just the $x$ direction (notice the second eigenvalue is 0).

Great, we got our principal components. Remember when I said we could fit this data as a linear combination of these two? To solve for the coefficients of these two vectors that fit the data, we simply have to use the Normal Equation again:

In [21]:
X = vec.copy()
B = np.dot((np.linalg.inv(np.dot(X.T,X))), np.dot(X.T,a.T))

fit = np.dot(vec, B)

plt.plot(x,y ,'ko', label='data', markersize=3)
plt.plot(fit[0], fit[1],'rx', label='PCA fit', markersize=2)
plt.legend()
plt.grid()
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Nifty. A couple of things to note - the number of principal componenents a matrix has is equal to its rank. Our matrix of $x$ and $y$ measurements was of rank 2, so there were 2 principal components. A data matrix can be fitted $\textit{exactly}$ with all of its principal components, including the noise. Of course, usually we don't want to overfit to noise, so in higher dimensional cases as we see below, you only typically need the first few components.

What we basically just did was rotate our initial vectors to match the data. Let's do another example where $x$ and $y$ are not so obviously correlated and the rotation is not so obvious:

In [22]:
phi = np.linspace(0, 2*np.pi, 501)

xx = np.cos(phi)+ np.random.normal(0, 0.1, 501)
yy = np.sin(phi)+ np.random.normal(0, 0.1, 501)

plt.plot(xx,yy,'ko', markersize=2)
plt.xlabel('x')
plt.ylabel('y')
Out[22]:
Text(0, 0.5, 'y')

While a circle seems meaningful, we should get in the habit of thinking of the variables $x$ and $y$ as columns in a data matrix $\mathbf{A}$:

In [23]:
A = np.column_stack([xx, yy])

fig, ax = plt.subplots()
im = ax.imshow(A, aspect='auto')
plt.colorbar(im)

ax.set_xticks([0, 1])
ax.set_xticklabels(['x','y'])
Out[23]:
[Text(0, 0, 'x'), Text(1, 0, 'y')]

We can do the same as above to solve for the principal components:

In [24]:
meanx = xx.sum() / N
varx = np.sum(( xx - meanx)**2) / N

meany = yy.sum() / N
vary = np.sum( (yy - meany)**2) / N

cov = np.sum( (xx - meanx) * (yy - meany) ) / N

K = np.array([[varx, cov],
             [cov, vary]])

val, vec = np.linalg.eig(K)

X = vec.copy()

B = np.dot((np.linalg.inv(np.dot(X.T,X))), np.dot(X.T,A.T))

fit = np.dot(vec, B).T


print(X.shape, B.shape)

fig, ax= plt.subplots()
plt.imshow(A, aspect='auto')
plt.title('Data')
plt.colorbar()
ax.set_xticks([0, 1])
ax.set_xticklabels(['x','y'])
plt.show()

fig, ax= plt.subplots()
plt.imshow(np.outer(B[0], X[0]), aspect='auto')
plt.title('component 1')
plt.colorbar()
ax.set_xticks([0, 1])
ax.set_xticklabels(['x','y'])
plt.show()

fig, ax= plt.subplots()
plt.imshow(np.outer(B[1], X[1]), aspect='auto')
plt.title('component 2')
plt.colorbar()
ax.set_xticks([0, 1])
ax.set_xticklabels(['x','y'])
plt.show()


fig, ax = plt.subplots()
plt.imshow(fit, aspect='auto')
plt.title('fit')
plt.colorbar()
ax.set_xticks([0, 1])
ax.set_xticklabels(['x','y'])
plt.show()

(2, 2) (2, 501)

Now let's apply this to our exoplanet data:

In [25]:
import pickle
wl_data, data, variance =pickle.load(open('Data/toi1431/20250701/cube.pic','rb'))
am = pickle.load(open('Data/toi1431/20250701/airmass.pic','rb'))
In [28]:
order = 25 #picking one at random

A = data[order].copy()
var = variance[order].copy()

n, m = A.shape #n exposures x m wavelength channels
wl = wl_data[order].copy()

fig, (ax1, ax2) = plt.subplots(nrows=2)

ax1.imshow(A, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')

#plot all spectra at once
ax2.plot(wl, A.T, color='r', alpha=0.05);
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())
Out[28]:
(1.7269065763762805, 1.7482534632516846)

I've organized my data into a 2D matrix of shape $n$ exposures/spectra $\times \ m$ wavelength channels. In the 2D image above, you'll notice that the counts are larger in the center of the spectrum. This is due to the instrument throughput and called the Blaze function. The dark, vertical stripes are absorption line from Earth's atmosphere (telluric features). Some of the weaker dark stripes might be stellar lines, although this particular exoplanet host star is a rapidly rotation A star, so there should be very few lines in this data. You also can't see it by eye yet, but on top of all this is some random noise. Of course, we want to get rid of all of these things so we can search for the faint exoplanet signal.

Even though each column is associated with a wavelength, and that wavelength has a physical meaning, think of this data matrix $\mathbf{A}$ as a collection of $m$ different random variables. Just like in the examples before, where the columns of the matrix were $\vec{x}$ and $\vec{y}$.

We want to find some $\hat{\mathbf{A}} = \hat{\beta} \mathbf{X} \approx \mathbf{A} + \epsilon$. How do we create the design matrix $\mathbf{X}$, and how do we choose the coefficients $\hat{\beta}$? There are a few ways we can go about this.

First, let's see if we can find the eigenvectors of the data covariance matrix like we have been doing:

In [29]:
# #first standardize to calculate the covariance matrix #
mean_in_time = A.sum(axis=0) / n
Astd = A - mean_in_time[None]
var_in_time = np.sum(Astd**2, axis=0)/ n
std_in_time = np.sqrt(var_in_time)
Astd /= std_in_time[None]

# cov = Astd.dot(Astd.T)
cov = Astd.T.dot(Astd)

val, vec = np.linalg.eig(cov) # eigendecomposition #

print(val)
print(vec)

[ 1.31291988e+05+0.j  1.10750965e+03+0.j  1.29895055e+02+0.j ...
  8.21739964e-14+0.j -1.72409823e-14+0.j  2.92589306e-15+0.j]
[[ 0.02268723+0.j -0.0030549 +0.j -0.04598958+0.j ...  0.0035415 +0.j
   0.02029295+0.j -0.02468113+0.j]
 [ 0.02312151+0.j  0.00476881+0.j -0.05801804+0.j ... -0.0015085 +0.j
   0.00218626+0.j -0.00252655+0.j]
 [ 0.0230276 +0.j -0.0042606 +0.j -0.06715401+0.j ... -0.00014659+0.j
  -0.00999681+0.j  0.01108111+0.j]
 ...
 [ 0.02332646+0.j  0.00983289+0.j  0.00630541+0.j ... -0.01338495+0.j
  -0.00693521+0.j  0.01959171+0.j]
 [ 0.02332568+0.j  0.01443843+0.j  0.00385698+0.j ...  0.02320953+0.j
  -0.0398242 +0.j  0.02735331+0.j]
 [ 0.02329988+0.j  0.01679966+0.j  0.00359763+0.j ...  0.05935766+0.j
  -0.01351007+0.j  0.02965148+0.j]]

Some of the eigenvalues and vectors are complex!! That's awkward since our data are 100\% real. (aside - this does work if you treat each exposure as the variables, rather than the wavelengths; i.e., PCA in the wavelength domain, but in general we want to work in the time domain to remove the constant in time signals.)

An alternative to eigendecomposition is something called singular value decomposition (SVD), which decomposes a matrix like so:

$$ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^\intercal $$

SVD can be preferable to eigendecomposition for a few reasons: 1) it is generally more numerically stable, 2) the vectors are all real, and 3) $\mathbf{A}$ does not have to be a square matrix, which means we can decompose our data matrix directly rather than dealing with the covariance matrix, which can be quite large ($m\times m$).

How to interpret the above expression? The matrix $\mathbf{U}$ is a $n \times n$ matrix whose columns are the $\textit{left singular vectors}$. There are $n$ left singular vectors because like eigendecomposition, the number of singular vectors is equal to the rank of the matrix. Each left singular vector also has $n$ elements, one for each of our exposures.

The matrix $\mathbf{\Sigma}$ is a diagonal $n \times n$ matrix. The diagonal elements are called the $\textit{singular values}$, which are related to and somewhat analogous to eigenvalues. Every singular vector has an associated singular value that is the effective weighting of how much variance that vector can explain, or how important it is.

Finally, the matrix $\mathbf{V}^\intercal$ is a $n \times m$ matrix of the $\textit{right singular vectors}$. There are $n$ right singular vectors associated with the corresponding left singular vector and singular value. Each right singular vector has $m$ elements, corresponding to one of our wavelength channels.

Let's see what all of that means in practice. Fortunately, there is a simple numpy function to do SVD:

In [30]:
LSV, S, RSV = np.linalg.svd(A, full_matrices=False)

print(LSV.shape)
print(S.shape) #since S is diagonal you only need a 1d array of the diagonal elements
print(RSV.shape)

(72, 72)
(72,)
(72, 1848)

Let's look at the first right singular vector:

In [31]:
plt.plot(RSV[0])
plt.show()

Hey, that kinda looks like our data! Let's see what happens when we multiply it by the first element of the corresponding left singular vector and the first singular value.

In [32]:
hmm = LSV[:,0][0] * S[0] * RSV[0]

plt.plot(wl, A.T, color='r',alpha=0.1);
plt.plot(wl, hmm, color='b')
plt.xlabel('Wavelength [micron]')
plt.ylabel('Counts [ADU]')
Out[32]:
Text(0, 0.5, 'Counts [ADU]')

That's a pretty good representation of the first spectrum in our data! If we plot the first left singular vector, we can see how the SVD wants to weight this right singular vector in time:

In [33]:
plt.plot(LSV[:,0])
plt.xlabel('exposure')
plt.show()

Coincidentally, this looks a lot like the mean count value in time. In fact, the first right singular vector also look a lot like the mean-in-time spectrum. It appears the SVD has identified that most of the variance in our data can be explained by just scaling the mean spectrum.

In [34]:
plt.plot(np.mean(A, axis=1))
plt.ylabel('Counts [ADU]')
plt.xlabel('Exposure')
plt.title('Meant Counts per Exposure')
plt.show()

plt.plot(wl,np.mean(A, axis=0))
plt.ylabel('Counts [ADU]')
plt.xlabel('Wavelength [micron]')
plt.title('Mean Counts per Wavelength')
plt.show()

Let's subtract this first singular vector from the data and see how well we've fit it:

In [35]:
first_vec = np.outer(LSV[:,0] * S[0], RSV[0]) 

fig, (ax1, ax2) = plt.subplots(nrows=2)

ax1.imshow(first_vec, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')
ax1.set_title('Raw Data')

ax2.plot(wl, first_vec.T, color='r', alpha=0.1);
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())

residuals = A - first_vec

fig, (ax1, ax2) = plt.subplots(nrows=2)
ax1.set_title('Data - SV1')
ax1.imshow(residuals, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')

ax2.plot(wl, residuals.T, color='r', alpha=0.1);
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())
Out[35]:
(1.7269065763762805, 1.7482534632516846)

Well that certainly fit most of the data, but we're left with these vertical stripes still. If you look closely, the wavelength channels of these stripes correspond to deep telluric lines. The depths of these lines vary with airmass and water vapor over the course of the night, so its not a surprise that we weren't able to fit them with just the median spectrum.

Maybe the second singular vector might help. Let's check:

In [36]:
plt.plot(wl,RSV[1])
plt.xlabel('Wavelength [micron]')
plt.show()

Seems to be mostly these troublesome telluric lines, so that's a good sign. What about the weights in time, i.e., the left singular vector?

In [37]:
plt.plot(LSV[:, 1])
plt.xlabel('Exposure')
plt.show()

Looks a lot like airmass!

In [38]:
plt.plot(LSV[:, 1] / np.std(LSV[:,1]), label='LSV 1')
plt.plot( -(am - np.mean(am)) / np.std(am), label='Airmass (standardized and negative')
plt.xlabel('Exposure')
plt.legend()
plt.show()

Let's subtract the second singular vector too:

In [39]:
scnd_vec = np.outer(LSV[:,1] * S[1], RSV[1]) 

fig, (ax1, ax2) = plt.subplots(nrows=2)

ax1.imshow(scnd_vec, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')
ax1.set_title('Data - SV1')

ax2.plot(wl, scnd_vec.T, color='r', alpha=0.1);
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())

residuals -= scnd_vec

fig, (ax1, ax2) = plt.subplots(nrows=2)
ax1.set_title('Data - (SV1 + SV2)')
ax1.imshow(residuals, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')

ax2.plot(wl, residuals.T, color='r', alpha=0.1);
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())
Out[39]:
(1.7269065763762805, 1.7482534632516846)

A lot better, but not perfect. It usually takes at least 4+ singular vectors to remove the strongest telluric features:

In [40]:
NPC = 5
weights = np.dot(LSV[:,:NPC], np.diag(S[:NPC]))
X = RSV[:NPC]

Ahat = weights @ X

residuals = A - Ahat

fig, (ax1, ax2) = plt.subplots(nrows=2)
ax1.set_title('Data - 5 SVs')
ax1.imshow(residuals, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')

ax2.plot(wl, residuals.T, color='r', alpha=0.1);
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())
Out[40]:
(1.7269065763762805, 1.7482534632516846)

But isn't that impressive?! With just a handful of vectors we are able to fit, in this case, 1848-dimensional data.

Before we move on, while you can do SVD on the raw data matrix, some people prefer to do so on the standardized matrix to make it "true" PCA. In this case, while the right singular vectors won't quite look like your data, the left singular vectors are true eigenvectors, which you can then use to fit for a suitable design matrix:

In [41]:
# standardize our matrix #
mean_in_time = A.sum(axis=0) / n
Astd = A - mean_in_time[None]
var_in_time = np.sum(Astd**2, axis=0)/ n
std_in_time = np.sqrt(var_in_time)
Astd /= std_in_time[None]

# do SVD #
LSV, S, RSV = np.linalg.svd(Astd, full_matrices=False)
NPC = 5 # pick a number #

B = np.column_stack([np.ones(n), LSV[:, :NPC]]) #add an intercept vector, bc the raw data is not mean subtracted 
X = np.linalg.inv(np.dot(B.T, B)).dot(np.dot(B.T, A)) #fit for vectors in wavelength

Ahat = B @ X

residuals = A - Ahat

fig, (ax1, ax2) = plt.subplots(nrows=2)

ax1.imshow(residuals, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')
ax1.set_title('Data - 5 PCs')

ax2.plot(wl, residuals.T, color='r', alpha=0.1);
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())
Out[41]:
(1.7269065763762805, 1.7482534632516846)

You can do the same thing by doing an eigendecomposition on the covariance between exposures (rather than between wavelengths):

In [42]:
cov = Astd.dot(Astd.T)

val, vec = np.linalg.eig(cov) # eigendecomposition #

B = np.column_stack([np.ones(n), vec[:, :NPC]])
X = np.linalg.inv(np.dot(B.T, B)).dot(np.dot(B.T, A))

Ahat = B @ X

residuals = A - Ahat

fig, (ax1, ax2) = plt.subplots(nrows=2)

ax1.imshow(residuals, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')

ax2.plot(wl, residuals.T, color='r', alpha=0.1);
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())
Out[42]:
(1.7269065763762805, 1.7482534632516846)

But like I said before, we generally want to do PCA/SVD in the time domain to remove signals that are stationary in time. Feel free to play with this on your own though.

ANYWAY - the fact that the second singular vector before corresponded to airmass variations was a "happy accident". There's no reason why airmass, and only airmass, has to be mapped onto one or only one of these vectors. But what if you "know" that airmass is an important variable and that the data variation should correlate with it? Or you at least want to start with a reasonable guess for the basis vectors of your data?

This is where SysRem comes in. Introduced by Tamuz+ 2005, the idea is that, if you have some good idea of how your data, $\mathbf{A}$, should vary with time, like airmass, in a vector $\vec{b}$, then you can solve for a new vector $\vec{x}$ that fits your data in each wavelength channel:

$$ \vec{x} = \begin{bmatrix} x_1, x_2, ..., x_m \end{bmatrix} \rightarrow x_j = \sum^n_{i=1} \frac{b_i a_{ij}/\sigma_{ij}^2}{b_i^2 / \sigma_{ij}^2} $$

If that expression looks familiar, that's because it's literally just the linear regression expression from way earlier, except weighted by the data variance $\sigma^2$ (as output by the spectral extraction).

Once you solve for $\vec{x}$, we subtract the 2D matrix $\vec{b} \otimes \vec{x}$ from $\mathbf{A}$, and then use $\vec{x}$ to solve for a new $\vec{b}$:

$$ b_i = \sum^m_{i=1} \frac{x_j a_{ij} / \sigma_{ij}^2}{x_j^2 \sigma_{ij}^2} $$

and iterate through this process several times.

In [43]:
niter = 5 #number of iterations

mean_in_time = A.sum(axis=0) / n

R = A - mean_in_time[None] #technically our first iteration is subtracting the mean ig

Barr = np.ones((n, niter+1)) #weights in time - the first is an intercept vector of 1's
Xarr = np.zeros((niter+1, m)) #vectors in wavelength
Xarr[0] = mean_in_time #set the first as the mean

B = am.reshape((n, 1)) #start with airmass

for i in range(niter):
    Barr[:, i+1] = B.reshape(n) #save it for later. why? its a secret ;)
    
    X = np.dot(R.T/var.T, B) / np.dot(1/var.T, B**2) #solve for x
    
    Xarr[i+1] = X.reshape(m) #save it
        
    R -= np.dot(X, B.T).T #subtract from existing residuals
    X= X.reshape(m)
    
    B = np.dot(R/var, X) / np.dot(1/var, X**2)# solve for new b
    B = B.reshape((n, 1))
    
    
Ahat =  Barr @ Xarr

fig, (ax1, ax2) = plt.subplots(nrows=2)

ax1.imshow(R, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')
ax1.set_title('Data - 5 SysRem Iterations')
ax2.plot(wl, residuals.T, color='r', alpha=0.1);
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())
Out[43]:
(1.7269065763762805, 1.7482534632516846)

Now, despite many people calling it as such, SysRem is techncially $\textit{not}$ PCA. Principal component vectors are orthogonal to each other. One way to check orthogonality is to take the dot product of two vectors. The dot product of two orthogonal vectors is 0.

In [44]:
print(Xarr[1].dot(Xarr[2]))

13763993.243948368

On the other hand, take our right singular vectors from SVD:

In [45]:
print(RSV[0].dot(RSV[1]))

-4.336808689942018e-17

What does this mean? In general, the features explained by a given singular vector do not correlate with features explained by another. But all of the variables captured by a given singular vector are correlated, so it's convenient to get all that in one fell swoop. This is not necessarily the case with SysRem - the different vectors can be correlated with each other, and a given vector might not capture only features that are correlated with each other. As in, the value of the vector in one wavelength channel doesn't necessarily depend on any other wavelength channel.

Take for example, the SysRem vector that is weighted by airmass with time:

In [46]:
plt.plot(wl, Xarr[1])
plt.xlabel('Wavelength [micron]')
plt.show()

Yes we have telluric lines as expected, by also a bit of the Blaze function, which physically is not correlated with the airmass at all.

In practice, this means that, perhaps for some wavelength channels, you need to do many iterations of SysRem, but for others, you only needed a couple, in which case you might be going a bit overkill in these channels and possibly corrupting the planet signal more than you need to. On the other hand, SysRem might capture some features that SVD missed because they didn't correlate with anything captured in the first few singular vectors. Your milage may vary.

At this point, you might be wondering how to choose the optimal number of components or iterations to remove. Unfortunately, there isn't really a good answer. Ideally, you want a number that successfully removes the telluric features, but does not eat away at your planet signal too much. One might be tempted to optimize your CCF signal strength, but this is likely to bias your inferences down the line, since you were optimizing to a specific model template on a specific velocity path for a specific noise instance. Personally, I think the best practice is to do your analysis several times with different number of components or iterations removed and make sure your answers match. Otherwise, it's just kinda based on vibes.

Speaking of accidentally removing some of the the planet signal, let's talk about model filtering.

As you might have gathered thus far, we don't really have a lot of control over what PCA, SVD, or SysRem map onto our new basis vectors/design matrices. We accept that to some degree, part of the planet signal has been captured by our detrending processes and removed, stretched, or otherwise corrupted. By not taking this into account, we could bias our inferences by overfitting our model to this disrupted planet signal. Fortunately, we have a pretty good idea of how the underlying planet signal was corrupted.

Let's read in a planet model spectrum and Doppler shift it in time:

In [47]:
from scipy.interpolate import splrep, splev
from scipy import constants
from scipy.ndimage import gaussian_filter1d

ph = pickle.load(open('Data/toi1431/20250701/ph.pic','rb')) #orbital phase
Kp = 188 #planet's radial velocity semi-amplitude
RV = Kp * np.sin(2 * np.pi * ph) #planet's radial velocity path

#read in a thermal emission model
wl_model, Fp, Fstar = np.loadtxt('papers/toi1431/figs/TOI1431_1DRC_FULL.txt').T

#I'm gonna smooth the star bc I don't feel like worrying about all the broadening nonsense rn...
Fstar = gaussian_filter1d(Fstar, 200)

#gaussian instrument profile 
xker = np.arange(41)-20
sigma = 5.5/(2.* np.sqrt(2.0*np.log(2.0))) 
yker = np.exp(-0.5 * (xker / sigma)**2.0)
yker /= yker.sum()

Fp = np.convolve(yker, Fp, mode='same')
Rp = 1
Rstar = 1

#spline interpolation coefficients
cs_p = splrep(wl_model, Fp*(Rp * 0.1)**2, s=0)
cs_s = splrep(wl_model, Fstar*Rstar**2, s=0)

We're going to initialize a $n \times m$ matrix $\mathbf{M}$ where each row is the Doppler shifted planet model spectrum. However, rather than saving as just $F_P$, we're going to save the planet model as $1+F_P/F_\star$. The idea is that, in theory our data is the combined signal of the planet and star's thermal emission multiplied by $\mathcal{T}$, a catch-all for the throughput/transmittance of our instrument and Earth's atmosphere, plus noise:

$$ \mathbf{A} = \mathcal{T} \times (F_\star + F_P) + \epsilon $$

Our detrending fit, $\hat{\mathbf{A}}$ shouldn't include the planet signal or noise, i.e., $\hat{\mathbf{A}} \approx \mathcal{T} \times F_\star$. In that case, multiplying $\mathbf{M}$ by $\hat{\mathbf{A}}$ will introduce our model planet spectrum in there already scaled by $\mathcal{T}$ and in the data units:

$$ \mathbf{M} \odot \hat{\mathbf{A}} = (1 + F_P/F_\star) \times \mathcal{T} \times F_\star = \mathcal{T} \times (F_\star + F_P) $$
In [48]:
M = np.zeros((n, m))
dl = RV * 1e3 / constants.c

fs_int = splev(wl, cs_s, der=0) #interpolate stellar spectrum at rest (or its RV) onto data wavelength grid

for j in range(n):
    wShift = wl * (1. - dl[j])
    fp = splev(wShift, cs_p, der=0) #interpolate planet model spectrum onto doppler shifted data wavelength grid
    M[j] = 1+(fp/fs_int)
    
fig, (ax1, ax2) = plt.subplots(nrows=2)

ax1.imshow(M, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')
ax1.set_title('Doppler Shifted Model Spectrum')
ax2.plot(wl, M.T, color='r', alpha=0.1, lw=0.1);
ax2.plot(wl, M[30], color='k', lw=1)
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('1 + Fp/Fs')
ax2.set_xlim(wl.min(), wl.max())
Out[48]:
(1.7269065763762805, 1.7482534632516846)

During our detrending process, we can save the final fit to the data matrix, $\hat{\mathbf{A}}$, as a noiseless representation of the data and inject the planet model into it. Then, with our original weights in time or eigenvectors, do the multilinear regression again and subtract the fit off once more. The idea is that we know that the detrending process is removing $\textit{something}$ from the planet signal and we know how that something varies in time.

For the SVD case:

In [49]:
NPC = 5
LSV, S, RSV = np.linalg.svd(A, full_matrices=False)
weights = np.dot(LSV[:,:NPC], np.diag(S[:NPC]))
X = RSV[:NPC]

Ahat = weights @ X #SVD fit

M_inject = M * Ahat #inject the planet model

X_new = np.linalg.inv(np.dot(weights.T, weights)).dot(np.dot(weights.T, M_inject)) #solve for new design matrix
Mhat = weights @ X_new 

M_filt = M_inject - Mhat #subtract off the fit

fig, (ax1, ax2) = plt.subplots(nrows=2)

ax1.imshow(M_filt, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')
ax1.set_title('SVD Filtered Model Spectrum')
ax2.plot(wl, M_filt.T, color='r', alpha=0.1, lw=0.1);
ax2.plot(wl, M_filt[30], color='k', lw=1)
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())
Out[49]:
(1.7269065763762805, 1.7482534632516846)

Funky. Let's take a closer look at how the model changed in a random exposure:

In [50]:
plt.plot(wl, M_filt[20]/np.std(M_filt[20]),label='Filtered Model (standardized)')
plt.plot(wl, (M[20] - np.mean(M[20])) /np.std(M[20]) , label='Initial Model (Standardized)' )
plt.xlim(1.735, 1.739) #zooming in on a small wavelength region
plt.legend()
plt.xlabel('Wavelength [micron]')
plt.show()

You'll notice that some of the planet lines are roughly the same, but others have completely different relative magnitudes or even signs!

We can do something similar after SysRem, which is why I saved each new $\vec{b}$ and $\vec{x}$ array as we iterated:

In [51]:
Ahat = Barr @ Xarr
M_inject = M * Ahat

Remember though that unlike SVD, SysRem fits to the variance-weighted data. Gibson+22 suggests adding in a whitening term:

$$ \hat{\mathbf{X}} = (B^\intercal \Lambda B)^{-1} B^\intercal \Lambda \mathbf{M} $$

where $\Lambda$ is a $n \times n$ diagonal matrix whose elements are the median inverse error per spectrum.

In [52]:
LL = np.diag( 1/np.median( np.sqrt(var), axis=1)  )

BB = LL.dot(Barr)
X_new = np.linalg.inv( np.dot(BB.T, BB)).dot( np.dot (BB.T, LL.dot(M_inject))  ) 

Mhat = Barr @ X_new

M_filt = M_inject - Mhat

fig, (ax1, ax2) = plt.subplots(nrows=2)

ax1.imshow(M_filt, aspect='auto', extent=[wl.min(), wl.max(), 0, n])
ax1.set_ylabel('Exposure #')
ax1.set_title('SysRem Filtered Model Spectrum')
ax2.plot(wl, M_filt.T, color='r', alpha=0.1, lw=0.1);
ax2.plot(wl, M_filt[30], color='k', lw=1)
ax2.set_xlabel('Wavelength [micron]')
ax2.set_ylabel('Counts [ADU]')
ax2.set_xlim(wl.min(), wl.max())
Out[52]:
(1.7269065763762805, 1.7482534632516846)
In [53]:
plt.plot(wl, M_filt[20]/np.std(M_filt[20]),label='Filtered Model (standardized)')
plt.plot(wl, (M[20] - np.mean(M[20])) /np.std(M[20]) , label='Initial Model (Standardized)' )
plt.xlim(1.735, 1.739)
plt.legend()
plt.xlabel('Wavelength [micron]')
plt.show()

Congratulations! You are now ready to perform cross-correlation or retrieval on your high resolution exoplanet data.